# 裁剪高分卫星的tif图像到矿区范围
import os

import numpy as np
from osgeo import gdal, gdal_array
from osgeo import gdal
gdal.UseExceptions()  # 启用异常处理
def is_coordinate_within_image(image_path, top_left_lon, top_left_lat, bottom_right_lon, bottom_right_lat):
    ds = gdal.Open(image_path)
    gt = ds.GetGeoTransform()
    x_min = gt[0]
    y_max = gt[3]
    x_max = x_min + gt[1] * ds.RasterXSize
    y_min = y_max + gt[5] * ds.RasterYSize

    if top_left_lon >= x_min and top_left_lat <= y_max and bottom_right_lon <= x_max and bottom_right_lat >= y_min:
        return True, x_min, x_max, y_min, y_max
    else:
        return False, x_min, x_max, y_min, y_max


def crop_image_by_coordinates(input_path, output_folder, top_left_lon, top_left_lat, bottom_right_lon, bottom_right_lat):
    in_ds = gdal.Open(input_path)

    # 获取输入图像的相关信息
    in_band1 = in_ds.GetRasterBand(1)
    in_band2 = in_ds.GetRasterBand(2)
    in_band3 = in_ds.GetRasterBand(3)

    # 获取原图的原点坐标信息
    ori_transform = in_ds.GetGeoTransform()
    top_left_x = ori_transform[0]  # 左上角x坐标
    top_left_y = ori_transform[3]  # 左上角y坐标
    w_e_pixel_resolution = ori_transform[1]  # 东西方向像素分辨率
    n_s_pixel_resolution = ori_transform[5]  # 南北方向像素分辨率

    # 计算给定经纬度坐标对应的像素位置
    top_left_offset_x = int((top_left_lon - top_left_x) / w_e_pixel_resolution)
    top_left_offset_y = int((top_left_lat - top_left_y) / n_s_pixel_resolution)
    bottom_right_offset_x = int((bottom_right_lon - top_left_x) / w_e_pixel_resolution)
    bottom_right_offset_y = int((bottom_right_lat - top_left_y) / n_s_pixel_resolution)

    # 确定切图的起始点像素位置和大小
    offset_x = min(top_left_offset_x, bottom_right_offset_x)
    offset_y = min(top_left_offset_y, bottom_right_offset_y)
    block_xsize = abs(bottom_right_offset_x - top_left_offset_x)
    block_ysize = abs(bottom_right_offset_y - top_left_offset_y)

    # 裁剪窗口超出图像边界的情况处理
    if offset_x < 0:
        block_xsize += offset_x
        offset_x = 0
    if offset_y < 0:
        block_ysize += offset_y
        offset_y = 0
    if offset_x + block_xsize > in_band1.XSize:
        block_xsize = in_band1.XSize - offset_x
    if offset_y + block_ysize > in_band1.YSize:
        block_ysize = in_band1.YSize - offset_y

    out_band1 = in_band1.ReadAsArray(offset_x, offset_y, block_xsize, block_ysize)
    out_band2 = in_band2.ReadAsArray(offset_x, offset_y, block_xsize, block_ysize)
    out_band3 = in_band3.ReadAsArray(offset_x, offset_y, block_xsize, block_ysize)

    output_filename = os.path.join(output_folder, os.path.splitext(os.path.basename(input_path))[0] + "_cut.tif")
    out_ds = gdal.GetDriverByName("GTiff").Create(output_filename, block_xsize, block_ysize, 3, in_band1.DataType)

    top_left_x1 = top_left_x + offset_x * w_e_pixel_resolution
    top_left_y1 = top_left_y + offset_y * n_s_pixel_resolution

    dst_transform = (top_left_x1, ori_transform[1], ori_transform[2], top_left_y1, ori_transform[4], ori_transform[5])
    out_ds.SetGeoTransform(dst_transform)

    out_ds.GetRasterBand(1).WriteArray(out_band1)
    out_ds.GetRasterBand(2).WriteArray(out_band2)
    out_ds.GetRasterBand(3).WriteArray(out_band3)

    out_ds.FlushCache()
    out_ds = None
    in_ds = None
def crop_images_in_folder(input_folder, output_folder, top_left_lon, top_left_lat, bottom_right_lon, bottom_right_lat):
    os.makedirs(output_folder, exist_ok=True)
    for filename in os.listdir(input_folder):
        input_path = os.path.join(input_folder, filename)
        if os.path.isfile(input_path) and (input_path.lower().endswith('.tif') or input_path.lower().endswith('.tiff')):
            within_bounds, x_min, x_max, y_min, y_max = is_coordinate_within_image(input_path, top_left_lon, top_left_lat, bottom_right_lon, bottom_right_lat)
            if within_bounds:
                crop_image_by_coordinates(input_path, output_folder, top_left_lon, top_left_lat, bottom_right_lon, bottom_right_lat)
            else:
                print(f" {filename} 裁剪区域超出该遥感图片，无法裁切，跳过该图")
                print("Image bounds: x_min={}, x_max={}, y_min={}, y_max={}".format(x_min, x_max, y_min, y_max))
# 示例用法
input_folder = 'E:\\坐标G6\\原始\\GF6_WFV_E105.8_N40.2_20231013_L1A1420364782'
output_folder = 'C:\\Users\\nxt\\Desktop\\数据集'

top_left_lon = 106.898  # 左上角经度
top_left_lat = 39.535  # 左上角纬度
bottom_right_lon = 106.963  # 右下角经度
bottom_right_lat = 39.445  # 右下角纬度

crop_images_in_folder(input_folder, output_folder, top_left_lon, top_left_lat, bottom_right_lon, bottom_right_lat)
# # 检查经纬度坐标是否在图像范围内，并获取最大范围
# # 检查经纬度范围是否在图像范围内，并获取最大范围
# within_bounds, x_min, x_max, y_min, y_max = is_coordinate_within_image(input_path, top_left_lon, top_left_lat, bottom_right_lon, bottom_right_lat)
# # 输出结果
# print("给定的经纬度范围：({}, {}) - ({}, {})".format(top_left_lon, top_left_lat, bottom_right_lon, bottom_right_lat))
# print("是否在图像范围内：", within_bounds)
# if within_bounds:
#     crop_image_by_coordinates(input_path, top_left_lon, top_left_lat, bottom_right_lon, bottom_right_lat)
# else:
#     print("给定的经纬度范围超出图像范围。图像的最大范围：")
#     print("x_min:", x_min)
#     print("x_max:", x_max)
#     print("y_min:", y_min)
#     print("y_max:", y_max)
#